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ABSTRACT 

On small angular scales, i.e. at high angular frequencies, beyond the damping tail of the primary 
power spectrum, the dominant contribution to the power spectrum of cosmic microwave back- 
ground (CMB) temperature fluctuations is the thermal Sunyaev-Zel'dovich (tSZ) effect. We investi- 
gate various important statistical properties of the Sunyaev-Zel'dovich maps, using well-motivated 
models for dark matter clustering to construct statistical descriptions of the tSZ effect to all orders 
enabling us to determine the entire probability distribution function (PDF). Any generic determin- 
istic biasing scheme can be incorporated in our analysis and the effects of projection, biasing and 
the underlying density distribution can be analyzed separately and transparently in this approach. 
We introduce the cumulant correlators as tools to analyze tSZ catalogs and relate them to corre- 
sponding statistical descriptors of the underlying density distribution. The statistics of hot spots in 
frequency-cleaned tSZ maps are also developed in a self-consistent way to an arbitrary order, to 
obtain results complementary to those found using the halo model. We also consider different beam 
sizes, to check the extent to which the PDF can be extracted from various observational configu- 
rations. The formalism is presented with two specific models for underlying matter clustering, the 
hierarchical ansatz and the lognormal distribution. We find both models to be in very good agree- 
ment with the simulation results, though the extension of hierarchical model has an edge over the 
lognormal model. In addition to testing against simulations made using semi-analytical techniques 
we have also used the maps made using Millennium Gas Simulations to prove that the PDF and 
bias can indeed be predicted with very high accuracy using these models. The presence of signifi- 
cant non-gravitational effects such as pre-heating, however, can not be modeled using an analytical 
approach which is based on the modeling of gravitational clustering alone. Our results indicate that 
the PDFs we construct are insensitive to the underlying cosmology and can thus provide a useful 
probe of non-gravitational processes e.g. pre-heating or feedback. 

Key words: : Cosmology- Sunyaev Zel'dovich Surveys - Methods: analytical, statistical, numeri- 
cal 



1 INTRODUCTION 

The inverse Compton scattering of C MB photons - known as the thermal Sunyaev-Zel'dovich effect (tSZ; ISunvaev & ZeldovichI i ll 9721 . Il980l) : 
iRephaelil ( fl99i) : iBirkinshavvi ( Il999l) ') - imprints a characteristic distortion in the Cosmic Microwave Background (CMB) spectrum that can be 
studied using surveys such as and the ongoing Planclfl satellite mission. The fluctuation of this distortion across the sky as probed 

by CMB observations can thus provide valuable clues to the fluctuations of the gas density and temperature. The up-scattering in frequency of 
CMB photons implies an increment in the spectrum at high frequencies with coixesponding decrement in the low frequency (Rayleigh- Jeans; RJ) 
regime, and a null around 217 GHz. This characteristic behaviour is a potential tool for the separation of tSZ from the other temperature anisotropy 
contributions. These techniques are extremely effective in subtraction of primary anisotropics due to its well understood (perfect black body) 

^ http://wmap.gsfc.nasa.gov/ 
^ http://www.rssd.esa.int/Planck 
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frequency dependence and near Gaussian statistical behaviour teouchet & Gispertlll999HDelabroullie^ Cardoso & Patanchodl2003l : iLeachlllOOSh . 
The tSZ effect is now routinely imaged in massive galaxy clusters where the temperature of the scattering medium can reach as high as lOkeV. 
This in effect produces a change in CMB temperature of order ImK at RJ wavelengths. 

Here we are interested in the general intergalactic medium (IGM) where the gas is expected to be at < IkeV in the sort of mild over-densities 
that lead to CMB contributions in the /iK range. In this work we primarily focus on the statistical study of wide-field CMB data where tSZ effects 
lead to anisotropics in the temperature distribution both due to resolved and unresolved galaxy clusters, keeping in mind that the thermal tSZ 
contribution is the dominant signal beyond the damping tail of the primary anisotropy power spectrum. We primarily focus on analytical modelling 
of the entire one and the joint two-point PDF of the tSZ effect. 

The tSZ power spectrum is known to be a sensitive probe of the amplitude of density fluctuations. Higher order statistics such as the skewness 
or bispectrum can provide independent estimates on the bias associated with baiyonic pressure, as well as providing further consistency ch ecks and 
cross - validation of lower order estimates. The modelling of lower order statistics of the tSZ effect done by various authors ( iSeliak.20od: Coorav 
20001 ; IZhang & Penll200ll : ICooravll2001bl ; lKomatsu & Seliakll200ll : IZhang & Setdl2007l) in the past has follow ed the halo model jCoorav & Seth 



2002 ) that relies on ingredients for the mass function based on the Press-Schechter ( Press & Schechteill 19741) formalism and radial profile given 

^ ays an important role in our understanding of the physic s involved 

jPersi et al.lll995Uda Silya et al.lll999l;lRefregier et al.ll2000l:ISel[ak et al.ll20od-|SDringel et alll200l|lwhite. Hemguist & Springelll2002l ; lLin et al.l 



l2004IZhang et alll2004 ICao. Liu & Fang|l2007l : lRoncarelli et alj|2007l ; iHallman et all200ll2009l : 1scott et alJl2012h .Though limited by their dy- 
namic range, some of these studies incorporates complication from additional radiative and hydrodynamical effects (i.e. "gastrophysics") such as 
radiative cooling, preheating and SN/AGN feedback to a certain extent which are otherwise difficult to incorporate in any analytical calculations. 

In parallel with the development o f this PS formalism, analytical modelling based on hierarchical form for the higher order correlation func- 
tions has also been studied ex tensively jSzapudi. Szalav & Boschailll992l ; ISzaDudi & Szalavlll993|;ISzapudi & Colombi|[T993 : ISzaDudi & Szalavl 
1 19971 : 1 Bemardeau et ail l2002h . We employ the particular form proposed by jBalian & Schaefferill989r) to model the tSZ statistics. This form 
has been studied extensively in the literature for modelling weak lensing as well for the statistics of collapsed objects and related astrophysi- 
cal phenomenon. The statistics o f collapsed objects and their contribution to the tSZ sky have been studied previously dVal ageas & Silk 199i; 
IValageas. Schaeffer & Sil3l2002h . These studies also probe X-ray luminosity from the same clusters. In this study we do not probe individual 
clusters or collapsed objects, but instead directly link the density field with corresponding SZ observables. Tho ugh the diffuse component of the 
tSZ e ffect is beyond WMAP detection threshold the situation may improve with future data sets such as P l anck ( Hansen et all 20051 ; [joudaki et al.l 
I2OIOI) or surveys such as Arc-minute Cosmology Bolometer AiTay Receiver (ACBAR); see lRunvan et alj |2003ln 

The paper is organized as follows. In Sj2]we provide the details of tSZ effect. We link the higher order multispectra of the SZ effect with 
the underlying mass distribution with the help of various biasing schemes in SjS] In ij4]we introduce the generic hierarchical ansatz and in Sj5]we 
introduce the specific formalism based on generating functions in the quasi-linear and highly nonlinear regime. In !j6]we show how the PDF and 
bias of the tSZ sky are related to that of underlying density PDF and bias. Various approximation schemes are discussed that can be used to simplify 
the PDF and bias. In S|7]we describe various simulations we have used in our study. In In Sj8]we present the results of tests agains simulations. 
Finally the ^is dedicated to the discussion of our result and future prospects. 



2 FORMALISM 

In this section we will provide necessary theoretical background for the computation of lower order moments of tSZ both for the one-point 
cumulants and the two-point cumulant correlators. These will be later used to construct the entire PDF and the bias of tSZ in the context of 
hierarchical clustering. We will be using the following form of the Robertson- Walker line element for the background geometry of the universe: 

rfs^ = -c^rff^ + a^(t)[rfr^ + d\(r){de'^ + sin ed<j>^)\. (1) 

Where we have denoted the comoving angular diameter distance by dA{r) and scale factor of the universe by a{t). dA{r) = K'^''^ sin{K+^'^r) 
for positive curvature, dA{r) = (— K)^'^''^ sinh((— K)"'"^''^r) for negative curvature and r for a flat universe. Here r is comoving distance or 
lookback time. For a present value of of Ho and f2o we have K = (fio + ^^a — l)Ho. The thermal Sunyaev-Zel'dovich (tSZ) temperature 
fluctuation ATgz (fi, i') = ST{Cl, v) /Tqms is given by the opacity weighted electron pressure: 

l:^Tsz{yi,v) = g^{x^)y{yi) = g^{x^) I dr ■jv^{Cl,r); 7re(x) = iSp<.(x)/(pe)- (2) 

Jo 

Here re is the Thomson optical depth; overdots represent derivatives with respect to r. Here y{Cl) is the map of the Compton y-parameter. 
The Thomson optical depth can be expressed in terms the Thompson cross-section ctt, is the Boltzman constant by the integral Te = 
c J ne{z)atdt. The free electron number density is represented by n^{z). The function gv{x) encodes the frequency dependence of the tSZ 



^ http://cosmology.berkeley.edu/group/swlh/acbar/ 
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Figure 1. Scatter plots of the baryonic pressure vs. the dark matter density contrast 1 + (5 = Pm/ (pm) are displayed for various simulations used in our study at a 
redshift 2 = 0. The left panels correspond to the simulations with no preheating (i.e. gravity only simulations but with adiabatic cooling), shown here as GO, where 
as the right panels correspond to the simulations with preheating, depicted as PC. More details about the simulations are presented in The smoothing scale for 
these plots is 0.5/i~^Mpc for the upper panels and 5/i~^Mpc for the lower panels. The scatter is lower for regions with high density contrast 1 + <5 > 100 where 
most of the tSZ signals originate. The simulations with preheating exhibits a less well defined correlation structure. We have shown 50,000 randomly sampled points 
from our simulations in each panel. 
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Table 1. The notations for various statistics of 5{r), ys(f!) and i/s(f^) are tabulated. The parameter y is constructed to have same statistics as 5 under certain 
simplifying approximation i.e. h{x) = h{x) and b{x) = b{x). The variance of S{r) is however different compared to that of j/s (fi). Also, notice that 5{r) is a 3D 
field whereas is a projected (or 2D) field. The normalised cumulants of y{fl) and <5(r) are identical, under certain approximation. Hence they are independent 
of the biasing. The variance of these fields are however different. 



anisotropies. It relates the temperature fluctuations at a frequency u with the Compton parameter y. Here the function gv{x) is defined as: 
gv{xi,) = coth (a::^/2) — 4 with = hv/ (kBTcMB) = i^/(56.84GHz). In the low frequency part of the spectrum gi,{xi,) = —2, 

for Xi, <^ 1, here x^, is the dimensionless frequency as defined above. We will primarily be working in the Fourier domain and will be using the 
following convention: 

7re(k;r) = J d'^x 7re(x; r) exp[— ik ■ x]. (3) 

The projected statistics that we will consider can be related to 3D statistics defined by the following expressions which specify the power spectrum 
P,r(fc; r), the bispectrum B,r(A:i, fe, fca; r) and the trispectrum r,r(fci, k2, ks, k^; r) in terms of the Fourier coefficients. In our notation we will 
separate the temporal dependence with a semicolon: 

(7re(ki;r)7re(k2;r))c = (27r)^5D(ki + k2)P^(fci; r); (7re(ki; r)7re(k2; r)7re(k3; r))c = (27r)^B^(fci, fe, fe; r)fo(ki + ka + ks) (4) 
(7re(ki; r)7re(k2;r)7re(k3; r)7re(k4; r))c = (27r)^r^(fci, fe, fcs, ^4; r)(5i5(ki + k2 + ka + k4). (5) 

The subscript c denotes the connected parts of a cumulant, i.e. those parts not related to the lower-order correlations. The multispectra of the 
underlying density field will develop a similar hierarchy. In the linear bias formalism, the multispectra of the field tt are directly linked to the 
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density contrast 5 — {pm — (pm) ) / (Pm) with pm being the h o mogeneous backgroun d density of the Universe. Using a bias factor 6,r('"), which 
depends on redshift, we can write dGoldberg & SDergellll999allbl; ICooravlbOOOl l2001bl) : 

P7r(fci;r-) = bl{r)Ps{ki;r); B,r(fci, ^2, fe; r^) = bl{r)Bs{ki, k2, k3;ri); r,r(fci, fc2, fe, fc4; r^) = bt{r)Ts{ki,k2,k3,k4;ri). (6) 

The redshift dependence of the bias is typically assumed to be of the following form: bTr{zs) ~ &,r(0)/(l + ^s) and 6^(0) — 
kBTe{Q)bs / {meC?). In our analysis we will show that it is possible to define a reduced y parameter whose statiscs will be insensitive to the 
details of the biasing model. The biasing scheme is motivated by a tight congelation of baryonic pressure and the density contarst at high density 
regions where most of the tSZ signals originate. 

We have studied the correlation between the fractional baryonic gas pressure Pg aB/(PBas) and the de nsity contrast pm/{/5?n) = 1 + 5 in 
Figure [T] The electronic pressure pe and the baryonic (gas) pressure Pgas are related (ICooravll200dll2001 alibi) p. = ||^Pgas; here X = 0.76 
is the primordial hydrogen abundance. The 3D electronic pressure fluctuation can be expressed as 1 + tt^ = Pe/(pe) = Pgas/(Pgas)- We 
find a correlation between these two variables i.e. 5 and ttg. The correlation is tighter for regions with higher density contrast 5, where most of 
the tSZ signals originate. We have considered two smothing scales. The top and bottom panels correspond to smoothing scales of .5ft^^Mpc and 
5/i~^Mpc respectively. We have displayed 50,000 randomly selected points from the simulations used for our study at a redshift of z = Q. The 
left panels correspond to the gravity only or GO simulation where as the right panels correspond to simulation with preheating and are denoted by 
PC. The correlation is more pronounced for the GO simulations. 



3 LOWER-ORDER STATISTICS OF THE THERMAL SZ EFFECT 



The statistics of the smoothed tSZ effect ys — {ys — (ys))/ (y) reflect those of the baryonic pressure fluctuations projected along the line of sight. 
Notice that in the denominator we have the average of unsmoothed y parameter. In our analysis we will consider a small patch of the sky where 
we can use the plane parallel (or small angle) approximation to replace spherical harmonics by Fourier modes. The three-dimensional electronic 
pressure fluctuations 7re(x) along the line of sight when projected onto the sky with the weight function a;sz(?') give the tSZ effect in a direction 
Q which we denote by y{Cl) (the smoothed y-maps will be denoted by ya{Cl) where subscript s will denote smoothed quantities): 



dr ujszir) ■k,{t,Ci); ujszir) = T,(r); y s{Cl) = dCl' Wci^ ~ Cl' ■,et)y{Cl') 



(7) 



Throughout, we will be using a Gaussian window Wg, or equivalently a Gaussian beam b; in the harmonic domain, specified by its full width at 
half maxima or FWHM — 6i,. Using the small ang le approximation we can compute the projected two-point correlation function in terms of the 
dark matter power spectrum Ps (k, r) l lKaisej|T992 



jfL exp(i e,2 ■ 1) bUot)P. 



I 



dA{r) 



ys = ys- (ys). 



(8) 



Here 612 is the angular separation projected onto the surface of the sky and we have also introduced 1 = dA{r)]s.± to denote the scaled projected 
wave vector. Using Limber's approximation, the variance and higher order moments of ys, smoothed using a Gaussian beam: bi (Ob) — exp[~l{l + 



with FWHM = 6b, can be written as: 



Jo dj^" '{r) 
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bi^iOb 



(27r) 



dA{r) 



dA{r) 



H + Ip 



(9) 



The higher-order moments of the smoothed temperature field relate {ys{db))c to the three-dimensional multi-spectra of the underlying pressure 
fluctuations B dHuHl 19991 : iMunshi & Colesll2000l) . We will use these results to show that it is possible to compute the complete probability 
distribution funct ion of ys from the underlyin g dark matter probability distribution function. Details of the analytical results presented here 
can be found in l lMunshi. Coles. Melotllll999j|bl: iMunshi. Melott. Coiesll2000l ; iMunshi & ColeslbOOd) . A similar analysis for the higher order 
cumulant correlators ( Szapudi & Szalavll997l : Mun shi. Melott. Coles'20001) of the smoothed SZ field relating (yf {di)yj({l2))c with multi-spectra 



of pressure fluctuation 

{rs{Cii)yu^2))c = 



can be expressed as (Munshi & Coles 200(1): 



.,P+9 



(r) 



(r) 



-dr 



X B 



(p+i) 



(27r)2 



bi^eb)--- 
fo(ii + 



dH. 



(2^)2 



{ei)bi^+, (eb) exp(z (li + ■ ■ ■ + Ip) ■ 6^2) 



(10) 



^dA{r) ' ' dA{r) 

We will use and extend these results in this paper to show that it is possible to compute the entire bias function b{> ys), i.e. the bias associated 
with those spots in the tSZ map where ys is above certain threshold, from the statistics of underlying over-dens e dark objects; this then acts as a 
generating function for the cumulant correlators. Details of the analytical results presented here can be found in dMunshi & Colesll200Cl) . 
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Dominant and Sub-dominant Contributions 




(c) (d) 
Dominant Sub-Dominant {fj= tf) 



Figure 2. A selection of diagrams tliat contribute to the cumulant correlators C21 and C31 are depicted. At the lowest order in non-Gaussianity, the cumulant 

coiTelator C'21 has two distinct contributions. The dominant contribution comes from the diagram (a). The other contribution from the diagram (b) adds negligible 

(2) —(2) (2) ~(2) 

contribution because of the multiplicative factor (^^ /^s )■ Here is the coiTelation function and is its average over the "cell" volume. We have suppressed 

the superscripts for brevity. The two diagrams denoted by (c) shows a representative dominant contribution to C31 . The diagrams denoted by (d) are the sub-dominant 

contributions to C31 . The upper panels of diagram (c) and (d) are of snake topology, where as the bottom panels correspond to the star topology. In our analysis we 

have ignored the sub-dominant term which are negligible when the cells are separated far apart. 



4 HIERARCHICAL ANSATZE : THE MINIMAL TREE MODEL AND ITS EXTENSION 

The spatial length scales corresponding to the small angular scales of relevance here are in the highly non-linear regime. Assuming a tree model 
for the matter correlation hierarchy in the highly non-linear regime for the density contrast 5, one can write the general form of the A'^th order 
correlation function S}^^ in terms of the two-point correlation function (}^^ jPeeblesI 1 980l ; IFtvII 1 984l Isemardeaijll 1 992h : 

(N-l) 

er(ri,...rN)= Qn,. n ^^'^(ri-rj). (H) 

a, N— trees labcllings edges 

It is i nteresting to note that a similar hierarchy develops in the quasi-linear regime in the limit of vanishing variance l lBemardeau & Schaeffeil 
I1992I) . However the hierarchical amplitudes Qn.q become shape-dependent functions in the quasi-linear regime. In the highly nonlinear regime 
there are some indications that these fu nctions become independent of shape, as suggested by studies of the lowest order parameter Q3 = Q using 
high resolution numerical simulations dScoccimarro & Friemanin"999l) . In Fourier space such an ansatz means that the hierarchy of multi-spectra 
can be written as sums of products of the matter power-spectrum: 

Bi(kl,k2,k3) = Q{P6{kl)Ps{k2) + Ps{k2)Ps{kA) + Ps{k3)Ps{kl))- (12) 

ra(ki,k2, ka, k4) = Ra Ps{ki)Ps{\ki + k2|)P«(lki + k2 + ka]) + cyc.perm. + Ri, Ps{ki)Ps{k2)Ps{k3) + cyc.perm. 

Different hierarchical mode ls differ in the way they predict the amplitudes of different tree topologies ( iMunshi. Coles. Melotj [l999allbl: 
iMunshi. Melott. Colesll2000l) . 

Working directly wi th cumulant generatin g function <j){y) (to be introduced later), modeling of the entire probability distribution function 
(PDF) was achieved by jColombi et"al]ll997l) . Remarkably, this was possible by a simple empirical modification of the existing quasilinear 
predictions and was dubbed Extended Perturbation Theory (EPT). The characteristic feature of this approach was to treat the local slope of 
the linear power spectrum as a free parameter in order to extend perturbative results to non-linear regime. EPT is fully parametrized by the 
non-linear variance and the sk ewness or equivalently Qs which is predicted by Hyper Extended Perturbation Theory (HEPT) developed by 
lacoccimarro &Friemanl(ll999h . 
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In recent years a new semi-analytical model for PDF has been proposed bv lValageas & Munshil JiooiFI . Its construction is very similar to 
EPT and works directly with 4>(y)- In addition to imposing the quasi-linear regime it also satisfies the rare-void limit. We will use this approach 
and show that the formalism is reasonably accurate in predicting the tSZ PDF. 

4.1 Cumulants 

Using this model we can express the one-point cumulants Sjv of ijs = ys — (ys) as: 

{yUOt)}c = (3Q3)C3[X|J = SsCslll] = Sz{yl{eb))l\ (13) 
{yt{Ob))c = {l2Ra + 4R,)C4ll] = SiCill'l] = S4y!{e,)}l (14) 

In general, for an arbitrary order cumulant we can write: 

i^smc = S,CA[XoJ-']; C, [[XeJ-'] = £^ ^^Jg^pej'^-^'dr; [I«J = J -^Ps (15) 

The projection effects on density c umulants were first derived in lBemardeau. Waerbeke. Mellieil(ll996h . In the context of weak lensing surveys the 
ea.l ll3t was derived bv lHuil ( fl999l) . Notice that Eq.dlSt relates the normalised one-point cumulants of ijs i.e. Sp to those of the underlying density 
contrast Sp. 

We will define a parameter ijs = ys/{y) which will further simplify the results. In particular we will find that, under certain simplifying 
assumptions, the statistics of y are identical to those of the underlying density contrast S. Notice that the definition of the parameter ijs involves 
the unsmoothed (ys) in the denominator. The normalised cumulants of ijs will be denoted as S = {y^)c/ {ys)c'~^\ which can be expressed as 

S = S{y)^^-'\ 

4.2 Cumulant Correlators 



The concept of cumula nt correlators was introduced by iBemardeaul ( 1996h . Later studies extended this result to observational studies including 



weak-lensing statistics jMunshi & Colesl200d : lMunshi & Jainl200d . l2001 ). In the present context, the family of cumulant correlators is important 



in modelling the bias associated with y maps. 

{fs{(li)ys{n2))c = 2QsC3ll0,Xe,4 = C2iC4XoJe,4 = C2i{fs)c{ys{^i)ysi^2))c, (16) 
{y'i{ni)ys{^2)}c = {3Ra + 6Rt)C4Xlle,4 = C3i{ys}c{ys{^i)ys{^2)}o. (17) 
In general, for an arbitrary order cumulant correlator we can write: 

(y?(ni)yj(n2))e = Cp,Cp+4[loJ+''-'lo,,] = CpM)i'^''~'''iysi^i)ysm)c, (18) 
where we have introduced the following notation: 

C,+,[[X«j(-+'-^)X9,J = £ [Io,4 = / ^Ps (^) exp(il ■ e..). (19) 

Notice that in the limiting case of 612 = we recover the limiting situation Cpg = Sp+q since we have [Xe^] = Xsia- The cumulant correlators 
for y denoted as Cpq and ij are related by the expression : Cpq = Cpq{yY^'^~^ . 

These lowest order statistics can be helpful in probing the pressure bias as a function of scale. It is however expected that the signal to noise 
will decrease with increasing order We will use these results to construct the entire PDF and the bias associated with the tSZ maps. This will be 
achieved using a generating function approach. 



5 THE GENERATING FUNCTION 

To go beyond order-by-order approach discussed so far we will use a formalism based on generating functions, which relies on the hierarchical 
scaling nature of the higher order correlation function. The scaling analysis deals directly with the generating functions that encode the information 
regarding the correlation hierarchy. The knowledge of these generating function are next useful in constructing the one- and two-point PDFs. 

In scaling analysis of the probability distribution function (PDF) the Void Probability distribution function (VPF) plays a most fundamental 

* http://ipht.cea.fr/Pisp/patrick.valageas/codepdf_en.php 
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Figure 3. The PDF ) is plotted as a function of {js . Vaiious curves con'espond to different beam size (FWHM) as indicated. Two different analytical models are 
shown; the lognormal (right-panel) and the hierarchical ansatz (left-panel). The solid lines correspond to the results from numerical simulations. The dashed lines in 
each panel represent the analytical results. 



role, because it can be related to the generating function (f>{z) of the cumulants or, if preferred, the Sp parameters l lWhitell979l : lBalian & Schaeffeij 
Il989l) : 

P„(0)=exp(-^), (20) 

— (2) 

where Py (0) is the probability of having no "particles" in a cell of of volume v, N is the average occupancy of these "particles" and Nc — N^g 

—(2) (2) 

and is the volume average of the two-point correlation function Q . Strictly speaking the above expression neglects any residual (subleading) 
Q dependence of the Sp parameters. The VPF is meaningful only for discrete distribution of particles and can not be defined for smooth density 
fields such as S or j/s(f2)- However the scaling function (j){z) defined above are very much useful even for continuous distributions where they 
can be used as a generating function of one-point cumulants or Sp parameters: (j)(z) = X^^i SpZ^ /p\. The function 4'{z) satisfies the constraint 
Si = S2 = 1 necessary for proper normalization of the PDF. The other generating function which plays a very important role in such analysis is 
the generating function for vertex amplitudes Vn, associated with nodes appearing in a "tree" representation of higher order correlation hierarchy 
{Q'A = V2,Ra ~ V2 and Rb = V'a)- In practice it is possible to work with a perturbative expansion of the vertex generating function Q{r). In terms 
of the vertices it is defined as: Q{t) = '^'^^Q{—T)^Vn/n\. However in the highly nonlinear regime a closed form is used. A more specific model 

for Q{r), which is useful to make more specific predictions (Bemardeau & Schaeffer 1992) is given by Q{t) = ^1 + t / K.a 

^{z) = zg{T)-\zT^g{r)- r = -z-^g{r). (21) 

—(2) ^ 

The range of S for which the power law regime is valid depends on the value of . For sm aller values of ^2 the power law regime is less 



pronounced. A more detailed discussion of these issues can b e found in Munshi et alj 1 1999|). T he links to the gravitational dynamics in the 



quasilinear regime, for various approximations are discussed in lMunshi. Sahni. Starobinskvl l ll994l) . The second equation of Eq. OTT l defines the 
variable z in terms of r; which can then be used to construct the function (piz) for a given model of C/(t). However a more detailed analysis is 
needed to include the effect of correlation between two or more correlated volume elements which will provide information about bias, cumulants 
and cumulant coiTelators of these collapsed object (as opposed to the cumulants and cumulant correlators of the whole map, e.g. iMunshi & JainI 
(I2OOQ, 2001)). We w ill only quote results useful for measurement of bias; detailed d erivations of related results including related error analysis c an 
be found elsewhere teemardeau & Schaeffejl 19991; iMunshi. Coles. Melotllll999j|bl ; IColes. Melott. Munshill 19991 ; iMunshi. Melott. Colesll2000l) 

Notice that t{z) (also denoted by I3{z) in the literature) plays the role of a generating function for factorized cumulant correlators Cpi 
{Cpq — CpiCqi): t{z) = y^^i Cpi/plz^. The PDF p(S ) and bias b(S) can be I' elated to their generating functions VPF 4i{z) and t{z) 
respectively by following equations ( iBahan & Schaeffeilll989l ; lBemardeau & Schaeffejl 99211 1999h : 

... r°° dz \ {l + S)z-4>{z) ^ r°° dz r(l + 5)£-0(£)] 
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It is clear that the function 4>{z) completely determines the behaviour of the PDF p{5) for all values of 5. The PDF too can be expressed with the 
help of a scaling functions: 



p{5) 



h{x) 



(1 + 5) 



The scaling function h{x) is related to (l>{z) through an inverse Laplace transform: 

hix) = — - — : / dzeyipixz^itjiz); h(x)b(x) 



(2) 



dz e-xp{xz)j3[z) 



(23) 



(24) 



The quantities and correspond to the density contrast 5 and $(2) and j3{z) will denote corresponding quantities for y^. The functions 
b{x) is simply a scaled version of 6(5): 



h{x) = 6 



1 + 5 

f{2) 
Si 



(25) 



We will next use these expressions to analyse the PDF and bias of y{Cl) maps. 



6 THE PDF AND THE BIAS OF BEAM SMOOTHED SZ EFFECT 



For computing the probability distribution function of the smoothed scaled tSZ field we will use the variable i/s that will make the analysis simpler: 
jjaidb) ~ ys{Ob) — {yaidb)), we will begin by constructing its associated cumulant generating function "l>(z): 



{y^{0b))c z" 



<i'(z) 



11 = 52 = 1. Now usin 



(26) 



Notice that the construction already satisfies the constraint equation 5*1 = 52 = 1. Now using the expressions for the higher moments of the SZ 
in terms of the matter power spectrum (see Eg. (115) ) gives: 



dA(r)2(f-i)(r) {yl{eb)}''J'-^'> 



We can now use the definition of (f>{z) for the matter cumulants to express ^(z): 



(27) 



(28) 



Note that we have used the fully non-linear generating function (j) for the cumulants, though we will use it to construct a generating function in 
the quasi-linear regime. Next we relate these results to the statistics of previously defined quantity ys. For the reduced tSZ field ys, the cumulant 
generating function ^{z) is given by, 



$(^) 



5^%! {y{e>,))Jo '^'Hyieb}). le. 



]cl>[b^{r){y{eb))- 



wsz{r) 



(29) 



The scaling_function h(x) for y associated with the PDF p(y) can now be related with the matter scaling function h{x) using the following 
definition ffialian & Schaeffer..l989l) : 



dz 

h{x) = ~ I - — :exp{xz)^{z) 



which takes the following form, using definitions corresponding to Eq, 



^ ^ Jo (ym 



{y{Ob)) Xe,wsz(r)6^(r) 



4(0 {ys{0b))c 



l0^wszir)b^{r) {y{9b)) 



(30) 



(31) 



^ Eq. )2U and Eq. )22t are cruci al to any approach that use s the hiera rchical ansatz for modeling of astroph ysical phenomena including the results presented here. 
These equations were derived bv lfiahan & Schaeffed l ll989h as well as lBemardeau & Schaeffel jl992lll999l) . 
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While the expressions derived above are exact, and are derived for the most general case using only the small angle approximation, they can be 
simplified considerably using further approximations. In the following we will assume that the contribution to the r integrals can be replaced by 
an average value coming from the maximum of LJsz(r), i.e. Vc (0 < Vc < rs). So we replace J f(r)dr by l/2/(rc)Ar where Ar is the interval 
of integration, and /(r) is the function of comoving radial distance r under consideration. Similarly we replace the wsz(r) dependence in the 1 
integrals by ajsz(»'c). Under these approximations we can write: 

^{z) = <i>{z)- h{x)^h{x). (32) 

Thus we find that the statistics of the underlying field 1+S and the statistics of the reduced field 5y are exactly the same under such an approximation 
(the approximate functions $(2) and h{x) do satisfy the proper normalization constraints). Although it is possible to integrate the exact expressions 
of the scaling functions, there is some uncertainty involved in the actual determination of these functions and associated parameters that describes 
it from N-body simulations (e.g. see Munshi et al. 1999, Valageas et al. 1999 and Colombi et al. 1996 for a detailed description of the effect of the 
finite volume correction involved in their estimation). 

In the following, we use $(2) as derived above to compute p{5y) with the help of Eq.ll22ll. In addition to the generating function approach 
we have used the lognormal distribution as a model for the underlying statistics (see Appendix-A for a detailed discussion on the lognormal 
distribution). 



6.1 The bias associated with the tSZ sliy 

The bias for the underlying density field S is defined using the joint two-point PDF p{Si, S2) for density contrasts 5i and S2 measured at two 
different points separated by a fixed distance. Given the two-point correlation function ^12 that characterizes the correlation hierarchy for this 
scale, the 2PDFp((5i, 82) can be expressed in terms of the one-point PDFs and the bias functions b{S) as follows: 

piSi,S2) = p{Sl)p(S2) [1 + fe(5i)Ci2fe(52)] . (33) 

The bias function b{5) is also useful in describing the bias associated with the over-dense regions. Cumulant correlators are the lower-order 
connected moments of the 2PDFp((5i, 52). 

To compute the bias associated with the peaks in the SZ field we have to first develop an analytic expression for the generating field l3{zi , Z2) 
for the SZ field ysifit) = ys{db) — {yiPb))- For that we will use the usual definition for the two-point cumulant correlator Cpq for the field: 

" {yl{9t)y+''-^{ys{^i)ysm)c' 

for a complete treatment of two-point statistical properties of smoothed fields see Munshi & Coles (1999b). We will show that, as is the case with 
its density field counterpart, the two-point generating function for the field ys can also be expressed (under certain simplifying assumptions) as a 
product of two one-point generating functions, /^^^(z), which can then be directly related to the bias associated with "hot-spots"in y-maps: 

Z2) = J2 T7^i ^1 - E E = ^(^i)/3(2/2) = r{zr)r{z2). (35) 

p,q p q 

The tree-structure of the correlation hierarchy that we have assumed is crucial to achieve the factorization derived above. It is also clear that the 
factorization of generating function actually depends on the factorization property of the cumulant correlators i.e. Cpg = CpiCgi. Note that such 
a factorization is possible when the correlation of two patches in the directions Cli and CI2 (j/s(f2i)j/s(n2))c is smaller compared to the variance 
{yt{Oh))c for the smoothed patches. The generating function (i{zi, Z2) for Cpq is constructed as follows: 

We will now use the integral expression for the cumulant correlators (Munshi & Coles 1999a) in order to express the generating function which, 
in turn, uses the hierarchical ansatz and the far-field approximation as explained above. Using Eg. (1191) we can write: 

/3(zi,^2) = y^ ^-1 ^2 ~ drdl{r)b^+''{r) [IsJ+^-'lg,,; (37) 

'^^ ^ p\q\ {yi{Ot)rc-' {yimn-' ^12 Jo ^' ^ ^dAw^pdAW^^' 

62 = {ys{(li)ys{^2)). (38) 

It is possible to further simplify the above expression by separating the summation over dummy variables p and q, which will be useful to establish 
the factorization property of two-point generating function for bias /3(2i, 22). We can now decompose the double sum over the two indices into 
two separate sums over individual indices. The above expression is quite general and depends only on the small angle approximation and the 
large separation approximation and is valid for any given specific model for the generating function G{t). However it is easy to notice that the 
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Figure 4. The cumulative bias associated with 6(> j/s) is plotted as a function of ijs- Various curves correspond to different beam size (FWHM) as indicated. Two 
different analytical models are shown; the lognormal (right-panel) and the hierarchical ansatz (left-panel), the solid lines correspond to the results from numerical 
simulations. The dashed lines in each panel represent the analytical results. 



projection effects as encoded in the line of sight integration do not allow us to write down the two-point generating function j3{z-i, Z2) simply as a 
product of two one-point generating functions /3(z) as was the case for the density field 1 + <5. As in the case of the derivation of the probability 
distribution function it simplifies matters if we use the reduced smoothed tSZ field . The statistical properties of ys are very similar to that of the 
underlying 3D density field (under certain simplifying approximations) and are roughly independent of the background geometry and dynamics of 
the universe, 

Jo {y[^b))^ £,12 le, V {yi{db))c d\{r) ) Xe, V {yi{Qb))c d\{r) ) 

While the above expression is indeed very accurate and relates the generating function of the density field with that of the tSZ field, it is difficult 
to handle in practice. Also it is important to notice that the scaling functions such as hix) for the density probability distribution function and h(x) 
for the bias associated with over-dense objects are typically estimated from numerical simulations specially in the highly non-linear regime. Such 
estimations are plagued with several uncertainties such as finite size of the simulation box. It was noted in earlier studies that such uncertainties lead 
to only a rather approximate estimation of h(x). The estimation of the scaling function associated with the bias i.e. h(x) (here a; = (1 + S)l^p^ 
is even more complicated due to the fact that the two-point quantities such as the cumulant correlators and the bias are more affected by finite size 
of the catalogs. So it is not fruitful to actually integrate the exact integral expression we have derived above and we will replace all line of sight 
integrals with its approximate values. Following our construction of the one-point PDF or p(5) we will replace integrals such as ^^^^ f{r)dr with 
their approximate values J^" f{r)dr = l/2/(rc)Ar; is at an intermediate redshift along the line-of-sight, its exact value is not important as 
the final result will be independent of Vc- For more accurate results we can Taylor expand /(r) or integrate the exact expression: 



(40) 



(41) 



Use of these approximations gives us the leading order contributions to these integrals and we can check that to this order we recover the factoriza- 
tion property of the generating function i.e. /3(2:i, 22) = I3{zi)l3{z2) = /3(zi)/3(22) = T{zi)f{z2). So it is clear that at this level of approximation, 
due to the factorization property of the cumulant correlators, the bias function b{x) associated with the peaks in the field {js, beyond certain thresh- 
old, obeys a similar factorization property too, which is exactly same as its density field counterpart. Earlier studies have established such a 
correspondence between the weak lensing convergence and the underlying density field in the case of one- and two-point probability distribution 
function p{S) (Munshi & Jain 1999b), 

b{xi)h{xi)b{x2)h{x2) = b{xi)h{xi)b{x2)h{x2). (42) 
Where we have used the following relation between P{z) = t{z) and b{x), 

b{x)h{x) — / dz t{z) exp(xz). (43) 

In Eq.(l32l( we have already proved that h{x) — h{x); hence from Eq.(l42t we deduce that b{x) — b{x). This means that the bias associated with 
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Figure 5. A simulated map of y-sky used in our study. The maps are 10 degree on a side and were generated on a 1024 X 1024 grid. These simulations were made 
using semi-analytical methods developed by iSchulz & Whitel2003h . More details about the simulations can be found in <Whitell2003h . 

dy is identical to that of the underlying density constrast 5. This is one of the main result of this paper. For all practical puipose we found that 
the differential bias as defined above is lot more difficult to measure from numerical simulations as compared to its integral counterpart where we 
concentrate on the bias associated with peaks above certain threshold, 

1 r°° t(z) 

b(> x)h{> x) = / dz exp(a;z). (44) 

It is important to notice that although the bias b[x) associated with the tSZ field ijs and the underlying density field are the same, the variance 
associated with the density field is veiy high but the projection effects in the tSZ field brings down the variance to a value comparable to unity. This 
indicates that we can use the integral definition of bias to recover it from its generating function (see Eq.(|43) and Eg. ( 144b ). Now, writing down the 
full two point probability distribution function for two correlated spots in terms of the tSZ field ysiBb)'- 

p(yi,y2)dyidy2 = p(yi)p(?/2)(l + b(yi)5f' (ri, r2)b(j/2))dj/idj/2; iji = ys(Oi). (45) 

The results derived here for auto-correlation can be generalised to cross-correlation analysis; this extension will be presented elsewhere. It is 
important to realise that the bias 6(ys) is related to the bias of "hot-spots" in the tSZ-map and relates their distribution with the overall correlation 
structure of the tSZ-maps. The bias 6^ defined earlier on the other hand relates the 3D pressure distribution 7re(x) to the underlying density 
distribution S{x). 

7 SIMULATIONS 

Two main approaches are generally used to simulate tSZ maps: (i) semi-analytical methods; (ii) Direct numerical simulations. In the following we 
discuss maps made using both sets of techniques. 

7.1 Simulations generated using semi-analytical Methods 

The semi-analytical methods identifies clusters from purely collisionless or N-body simulations; simplifying assumptions regarding the distribu- 
tion, hydrodynamical and thermodynamic equilibrium properties of the baryons are then used to create a baryonic data cube from the underlying 
dark matter distribution. The resulting distribution of baryons is then eve ntual ly used to cr eate tSZ y maps. This method was pioneered by 
iKav. Liddle & Thomas! fiooH) and was later used bv lSchulz & White! ( 12003!) and !White! ( !2003!) . In our study we have made use of data described 
in (White 20(^ where more detailed discussion about these simulations as well as the process of generating the tSZ y-maps that we have used 
can be founclf]. The numerical simulations are rather expensive and semi-analytical methods can provide reasonably accurate results that can be 
used to understand the underlying physical processes which affect the statistics of tSZ process. The particular cosmology that we will adopt for 
numerical study is specified by the following parameter values : S^a = 0.741, h — 0.72, Q,t — 0.044, JIcdm = 0.215, ris — 0.964, as ~ 0.803. 

® http://mwhite.berkeley.edu/tSZ/ 
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Figure 6. Simulated 5° X 5° dimensionless scaled thermal Sunyaev-Zel'dovich maps logj^(,[?// (j/)] are depicted. The maps were generated using Virgo consortium's 
Millennium Gas Simulation. The left panel shows the resulting y maps. The middle panel coirespond to maps generated using low density regions. Only over dense 
regions with density 1 + 5 < 100 were considered. The right panel correspond to low temperature regions T < lO^K. These set of hydrodynamic simulations 
ignore pre-heating but takes into account adiabatic cooling. We will refer to them as GO or Gravity-Only simulations (see text for more details). 




Figure 7. Same as previous figure, but for simulations with pre-heating and cooling. These simulations will be referred to as PC simulations. 



7.2 Hydrodynamic Simulations 

Early simulations of the tSZ modeled gravitational collapse and adiabatic compression but initially ignored the effect of adiabatic cooling which 
affects the thermal state of the gas especially in the halos. It was soon realised, however, that the baryons are subjected to non-gravitational heating 
processes - such as the feedback of energy from supemovae or AGN. Over the years the field of numerical simulations has matured to such an 

extent that it is now possible to include these effects self-consistently and in a r easonably large simulation box. 

The simulated y-maps that we have used were generated by IScott et a l. (2012) using millennium gas simulations jHartlev et al.l |2008| ; 
IStarek. Rudd & Evrardll2009t lYoung et al.ll201ll: IShort et al.lbOlOl) . Which in turn were generated to provide hydrodynamic versions of the Virgo 
conso rtium's dark matter Millennium Simulations and were performed using publicly-available GADGET2 N-body/hydrodynamics code jSpringell 
l2005h . Two different versions of the simulations use same initial conditions and box-size. In the first run, the gas was modelled as ideal non-radiative 
fluid and was allowed to go adiabatic changes in regions of non-zero pressure gradient. The evolution was modelled using smooth particle hydro- 
dynamics (sph). An artificial viscosity too was used to convert bulk kinetic energy of the gas into its internal energy. This is essential to capture 
the physics of shock and thus generate quasi-hydrostatic equilibrium. These process ensures quasi-hydrostatic equilibrium inside vitalized halos. 
See text for more details of the hydrodynamic simulations used to generate these maps. These set of simulations will be refen' ed as Gravi ty Only 
(GO) simulations. Non radiative descriptions of inter-cluster gas do not reproduce the observed X-ray properties of the clusters ( IVoitl2005l) . So the 
next set of simulations that we use pre-heated gas at high redshift that can generate the required core entropy and capable of producing a steeper 
X-ray luminosity-temperature in agreement with obse rvations . The entropy level of these sec ond set of simulations were chosen to match the 
mean X-ray luminosity temperature relation at z = l lKaiser n] |1991| : lEdward&Henrvlll99ll) . These simulations also include radiative cooling 
and an entropy sink. We will refer to these simulations as PC. Cooling in these simulations do not play an important role as the cooling time 
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Gravity and Adabatic Cooling Simulations: PDF 




1+y 1+y 1+y 

Figure 8. The probability distribution func tion for y is being compared to theoretical predictions from lognormal distribution and an extension of hierarchical model 
as proposed by iValageas & Munshil2004l). The solid-lines repre sents the PDF computed from the simulation. The soHd and dashed lines coiTespond to the lognormal 
model and the hierarchical model of l lValageas & Munshill2004 . The smoothing angular scales correspond to Bq = 0.25" (left-panel) 80 = 1.25' (middle-panel) 
and 80 = 2.5' respectively. The numerical curves are averages of three individual realisations each. 



for the preheated gas is long compared to the Hubble time. The cosmological parameters of these simulations are JIcdm = 0.25, JIa ~ 0.75, 
S^i, = 0.045, h = 0.73 and erg = 0.9. 

The scaled logjg[j//{2/)] parameter distribution of a realisation is shown in Figure ^ for gravity only or GO simulations and Figure (|7] for 
simulations with pre-heating and cooling (PC). The left panels show contribution from all individual components. The middle panels represents 
contribution from over dense regions that satisfy the constraint 1 + S < 100. Finally the right panels correspond to the contribution to the y-map 
from gas which satisfy the constraint T < lO'^K. 

There is a very clear and obvious difference between the two sets of maps in that the GO maps have more substructure. The smoothness of 
the PC maps is due to the external thermal energy added to the gas by the pre-heating process. The mean y-parameter in the GO simulations is 
(y) — 2.3 X 10~® and in the PC simulations it is nearly four times higher (y) = 9.9 x 10~®. These values are consistent with COBE/FIRAS 
constraint (y) < 1.5 x 10^^. However, it is believed such a high level of preheatin g would definitely rernove some of the absorption f eatures seen 
in the Lyman-a spectrum observed towards quasars dTheuns. Mo. Schave 2001: Shang. Crotts. Haiinanll2007l : Isorgani & Vielll2009l) . Indeed the 
PC model studied here should be treated as an extreme test of the effect of a high pre-heating scenario. 

In terms of source contributions, the bulk of the y-signal comes from low redshift i.e. z < 2. However in case of PC simulations the opposite 
is true, where 80% of the signal originates from z < 3.5. The overdense regions such as the group or clusters are the sources of y-signal in the 
GO simulations which are primarily embedded in structures that collapsed at relatively lower redshift. In case of the PC simulations most of the 
signals comes from mildly overdense gas at high redshift. It's interesting to notice that the GO simulations do get contributions from the gas at 
high redshift 2 > 4. However such contributions are completely erased in case of the PC simulations. This is primarily due to the fact that radiative 
cooling erases most of the ionized gas at these redshifts. 



8 TESTS AGAINST NUMERICAL SIMULATIONS 

In this section we present the result of our comparison of theory against the simulations described in the previous section. We have checked our 
results against simulation that are based on N-body simulations where baiyon is added using a semi-analytical prescription. We have also compared 
our results against state of the art hydrodynamic simulations. 



8.1 Semi-analytic Simulations 

We have used simulated y-maps described in l lWhitell2003h to test our theoretical predictions. The simulations were generated on a 1024 x 1024 
grid and cover 10° x 10° patches on the surface of the sky. To compute the PDF and the bias we construct the reduced {js maps from the beam- 
smoothed j/s- A Gaussian beam with varying FWHM 6t = 30", 1', 5', 10' was used for this purpose. The binning of the data at each grid point 
was next performed to the histogram and finally the one-point PDF of j/s • In Figure (|3} we have presented the results of our calculations for the 
entire range of beam sizes and compared them against the results from the simulations. 

For computation of the bias b{ys) from the precomputed ya we found that the estimation of cumulative bias 6(> j/s) is much more stable 
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Gravity and Adiabatic Cooling Simulations: PDF 
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Figure 9. The probability distribution func tion for y is being compared to theoretical predictions from lognormal distribution and an extension of hierarchical model 
as proposed by iValageas & Munshil2004l). The solid-lines repre sents the PDF computed from the simulation. The soHd and dashed lines coiTespond to the lognormal 
model and the hierarchical model of l lValageas & Munshill2004 . The smoothing angular scales correspond to Bq = 0.25" (left-panel) 6o = 1-25' (middle-panel) 
and 00 = 2.5' respectively. The numerical curves are averages of three individual realisations each. 
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Gravity and Adiabatic Cooling Simulations : Bias 
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Figure 10. The bias fe(> j/s) plotted as a function of ya for various smoothing angular scales correspond to da = 0.25" (left-panel) Bq = 1.25' (middle-panel) and 
do = 2.5'. Theoretical predictions from hierarchical ansazt (short -dashed), lognormal (long-dashed) and simulations (solid-lines) are shown. 



then its differential counterpart b(ys)- We used the following expression for computation of the bias function: 

dyi 



b{> Vs) = 



/oc 
is 



dy2p{yi,y2) ^ 
[J^ dyip{yi)]'2 



(46) 



The one- and two-point PDFs as well as the two-point correlation function are all constructed at grid points. We found that the numerical compu- 
tation of b{> ys) to be much more stable than that of the differential bias b{ys). The computation of two-point PDFs were done using different 
separation angular scales 9i2. The theoretical bias is computed at the large separation limit. However, we found that this limit is reached very 
quickly, typically at an angular scale which is twice the FWHM i.e. 612 ~ 26o. We also notice that the estimation of the bias is independent of any 
assumption regarding the factorizability of 2PDF; i.e. Eq.J46|) do not depend on any such ad hoc assumptions. 

In Figure (|4j we have shown the integrated bias b{> ys) computed using two different analytical techniques as function of the threshold ys- 
In the right panel the various smoothing angular scales FWHM we have shown correspond to 9t — 30", 1', 5', 10'. The solid lines correspond to 
the numerical simulations. The dashed lines correspond to predictions from the lognormal model. This bias for higher {js regions is higher. The 
results for the perturbative calculations are shown only for di, = 10', 5'. For smaller beam the perturbative calculations break down. 

By construction, the statistics of the y field are insensitive to the background cosmology, in that they mimic the statistics of the underlying 
density contrast S. Hence the shape of its PDFs do not depend on the detailed modelling of the electron pressure or the related electron-pressure 
bias b,r{r). The variance of these distributions however do depend on this function. We have used the variance {j/s)c computed from the maps in 
our computation of the PDFs. 
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8.2 Hydrodynamical (SPH) Simulations 

In addition to the maps generated using semi-analytical methods we have also used maps generated by realistic state of the art hydrodynamic 
simulations. These maps are 5° x 5° in size and are constructed using a 1024 x 1024 grid. We have analysed two different sets of maps to test our 
analytical predictions which we describe below. 



8.2.1 Simulations with Adiabatic Cooling: 

The first set of maps are derived from th e GO simulations described previously. To these we have compared both the lognormal approximation and 
the predictions from analytical results o f IValageas & MunshilJlOOj) . We have tried three different angular scales 0(, = 1 5" , 1 . 25' , 2 . 5' respectively. 
We numerically evaluate the statistics of ys, using exactly the same technique described before. We find that the results from numerical simulations 
are reproduced extremely well in the realistic and up-to-date hydrodynamic simulations. In case of the GO simulations, where gravitational 
dynamics and adiabatic cooling are the main factors influencing structure formation ,the numerical PDF is reasonably reproduced by theoretical 
predictions. We have presented these results in Figure jS). As noted before, perturbative predictions start to break down when the variance at a 
given scale reaches unity. The theoretical predictions of Valageas & MunshU ( |2004|) remedies the situation and is useful for analytical prediction 



of PDF at an arbitrary non-linear scale. We find that analytical prediction for ijs is accurate down to very low values of PDF. It reinforces our 
conclusions drawn from analysis of maps generated using semi-analytical techniques. Though the maps from hydrodynamical simulations have 
less sky-coverage compared to the semi-analytical maps used before they have more realistic representation of the baryonic physics responsible 
for the tSZ effect. 

In addition to studying the maps, we have divided the entire contribution from various baryonic components, to check how our theoretical 
prescriptions compare with that from simulations for individual components. In this context, we notice that, thermodynamic states of baryons, 
as well as their clustering, at low to med ium redshift z < 5, has been studied, using both numerical as well as analytical techniques. In their 
studies, IValageas. S chaeffer & Silkl ( |2002|) has used the hierarchical ansatz, to study the phase-diagrams of cosmological baryons as function of 
redshift. The low temperature "cool" component of the intergalactic medium (IGM) represented by Lyman-Q forest typically satisfies the constrain 
lO^K < T < 10*K. The exact values of the lower and upper-limit depends somewhat on the redshift. The "cold" component of the IGM is very 
well characterized by a well-defined equation of state. The "warm" component of the IGM on the other hand is shock-heated to a temperature 
range of lO^K < T < lO^K due to the collapse of non-linear structure and can not be defined by a well defined equation of state. Though 
the "warm" component does follow a mean temperature-density relation, the scatter around this relation however is more significant than for the 
"cool" component. Both the "cool" and "warm" components originate outside the collapsed halos and typically reside in moderate overdensites 
1 + (5 < 100. Finally the remaining contribution comes from the hot baryonic component of the virialized high density halos with temperatures 
T > lO^K. 

In Figure ([9} we have compared the contributions from the medium over-density regions 1 + 5 < 100 (medium panel) which is caused by 
both "warm" and the "cool" component of the IGM. The right panel correspond to emission primarily from the "cool" component. The variance 
corresponding to these individual components are computed individually and used as an input in the computation of their PDF. The results are 
presented for the particular case of smoothing 6s — 15" although the agreement is equally good for other smoothing angular scale. The sharp 
drops seen in the PDFs are due to the finite size of the catalog. The lowest probability that we can compute using a 1024 x 1024 grid is roughly 
10~®. For larger smoothing angular scales it is slightly less. 

The computation of integrated bias h{> ys) was carried out using techniques discussed in the previous section. We find that the analytical 
predictions for bias to match very accurately with the one recovered from numerical simulations. These results are presented in Figure The bias 
were estimated for the same angular scale. Indeed, in almost all cases the predictions for lognormal as well as hierarchical model are very close. 
An extension of lValageas & Munshil ( |2004|) technique for the case of bias is possible, but has not been worked out. However, bias computed using 
hierarchical model and lognormal approximation seems to be reasonably accurate. The bias 6(> ys) being a two-point statistic is more sensitive 
to sample variance. 



cS.2.2 Simulations with Adiabatic Cooling and Pre-Heating: 

The second set of simulations that we have analysed includes pre-heating and is referred to as the PC simulations. The smoothness of the PC maps 
are reflected in their low variance. This is primarily due to the high level of the pre-heating that erases many substructures resulting in maps with 
less features. We include these simulations in our studies mainly to test the limitations of analytical predictions. The fundamental assumption in 
our analytical modeling is of gravity induced structure formation where baryons are considered as the biased tracers of underlying dark matter 
clustering. We find significant deviation of the numerical results from theoretical predictions in the presence of high level preheating at small 
angular scales Figure ( III) . These deviations are more pronounced at smaller angular scales. The PDFs become Gaussian at scales 9s ~ 10' or 
larger. The deviation at all-scales is less pronounced if we remove the collapsed objects and focus on the maps with overdensites 1 + (5 < 100. 
However the PDF of the y distribution from "cold" intergalactic gas is represented very accurately by our analytical results at all angular scales in 
the presence of pre-heating Figure (|9). 
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Figure 11. The PDF p(y ) of y is compared to theoretical predictions from lognormal distribution and an extension of hierarchical model as proposed by 
jValageas & MunshillOM) . From left to right panels correspond to total contiibution, contribution from regions of low overdensity 1 + 5 < 100 as well as 
contribution only from low temperature regions T < 10''' K of the simulations. In each panel we show two smoothing scales 8s = 1.25' as well as Bs = 5' 
respectively. 



9 CONCLUSIONS 

We have studied the prospects for extracting and using the non-Gaussian statistical signatures from tSZ maps. The tSZ effect is associated with the 
hot gas in large-scale structure that is probed by multi-frequency CMB experiments. When compared to the CMB temperature anisotropics, the 
tSZ effect has a distinct spectral dependence with a null at a frequency of 217 GHz. This distinct spectral signature of SZ effect means it can be 
effectively separated from the primary CMB contributions. This will provide an unique opportunity to probe tSZ effect using data from ongoing 
surveys. 

The statistical analysis of frequency-separated tSZ maps has so far been mainly focused on lower order statistics or topological descriptors 
jMunshi et al.l2012l) . Non-Gaussianity in the tSZ signal is an additional information that is useful in constraining the large scale pressure fluctuation 
associated with the tSZ effect. This signature can be useful for constraining non-gravitational effects such as pre-heating or other forms of energy 
injection in the form feedback from active galactic nuclei (AGN) or super novae (SN). 

The tSZ effect traces the pressure fluctuations associated with the large scale distribution of the baryonic gas. The pressure fluctuations due 
to the virialized dark matter halos can be modeled by assuming them to be in hydrostatic equilibrium with the dark matter distribution in the halo. 
Such a halo model des cription has been extensively used in understanding the statistical properties of tSZ effect as well as other CMB secondaries 
l lCoorav & Seti]|2002l) . 

In addition to the halo model prescription, a redshift-dependent linear biasing sche me that relies on a perturbativ e description of dark matter 
clustering has also been in use to model certain aspects of tSZ e ffect. It was introduced by Goldberg & Spergell(ll999allbl) . Combining such a model 
with prescriptions from hyper-extended perturbation theoiy of lScoccimarro & Frie man Vl999^. can be used to predict the lower-order moments; 



this approach has been tested successfully against numerical simulations bv lCoorav et al. (2000D. 

In the following we summarize the main conclu sions of this study. 

Use of pe rturbative approach and its extension bv lValageas & MunshJ j2004t) : We have used the approach developed bv lGoldberg & Spergell 
jl999allbh : ICoorav et al.l 1 I2OO0I) to construct the PDF and bias. We construct the entire cumulant-generating function and show that under certain 
simplifying assumptions it becomes independent of the details of the biasing scheme. The generating function adopted here was developed pri- 
marily for the construction of 3D and 2D (projected) density distribution that are studied using galaxy surveys; later it was used in construction of 
weak lensing PDF in small and large smoothing angular scales. In this paper we have shown that a similar technique can be adopted for the study 
of tSZ when it is expressed in terms of the 3D pressure fluctuations. We have compared these results against two sets of maps: (a) maps made from 
semi-analytical simulations which are of 10° x 10°in size (b) 5° x 5° maps using full SPH simulations (the millennium gas simulations from Virgo 
consortium). We have studied angular scales 10" < 61, < 10' for the semi-analytical maps and 2.5" < 61, < 5' for the smaller maps. On angular 
scales where the rms fluctuation in ys maps is lower than unity the perturbative series is valid and the analytical predictions are in very good agree- 
ment with numerical simulations; we find analytical results from perturbation theory to be accurate for angular scales larger than few arc-minutes 
(6b > 2'). As a result of our study we notice that even at comparatively large angular scales, 9t > 10', the PDF distribution of th e tSZ effect is 
highly non-Gaussian. Going beyond the perturbative approach we have used the formalism presented in IValageas & Munshil ( l2004h to extend the 
analytical predictions to all angular scales. We find an excellent agreement with theory and simulations for all possible angular scales studied by 
us for simulation where gravity and adiabatic cooling plays a dominant role. Use of lognormal approximation: In addition to the perturbative 
approach and its extension using hierarchical ansatz we have also used a model based on the lognormal distribution. The lognormal model is non- 
perturbative and has been used widely in the literature. Like the perturbative approach it has been used to model the results from galaxy surveys 
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T<10^K systems in Pre — Healing Simulations: PDF 




1 1.5 0.8 1 1.2 1.4 0.8 1 1.2 

1+^ 1+y 1+y 

Figure 12. The pdf for y, p{y) from numerical simulations (solid-lines) for "c old" T < 10^ K com ponents are compared against the theoretical predictions from 
lognormal distribution (dashed-lines). The hierarchical model as proposed by jValageas & Muns "hi| [2004!) produces nearly identical results and are not shown. The 
angular scales considered from left to right are 9s = 0.25", 9s = 1.25' and 9s = 2.5' respectively. 



faamiltodl 19851 : IColes & Jon ej'l99lVBouchet et all l 19931 : iKofman et al.lll994h . weak lensing obsei^ables jMunshillloOOl : IXaniva et alj|2002l) as 
well as Lyman-alpha statistics (Bi & Davidson 19 97f). The validity of the lognorm al model has been compared against numerical simulations as 
well as against peiturbative or hierarchical methods (lBemardeau & Kofman|[l993) . In the entire range of angular scales considered by us we find 
the lognormal model to be a very good approximation for the y parameter distribution. We have also used the bias computed from lognormal 
approximation and found it to be in reasonable agreement. 

Isolating the Effect of Background Cosmology: We have shown that the statistics of y can be described using analytical models of gravitational 
clustering alone and the Sn parameters that describe the PDF as well as the entire PDF is insensitive to the background cosmology - the lognormal 
model e.g. do not have any cosmology built into it. The power spectmm of y on the other hand is sensitive to the amplitude of the density 
fluctuations, erg and other cosmological parameters. Using higher order statistics , such as the skewnes s, previous authors argued the possibility 
of separation of the pressure bias from the amplitude of the density fluctuations jffill & Sherwinll2012h . The dimensionless parameter y that we 
have introduced here is insensitive to background cosmology and in this sense our approach achieves separation of background cosmology and the 
effects of gravitation to all order. 

Separation of gravitational and non-Gravitational Aspects: Though the PDF constructed using the semi-analytical approaches do agree with 
numerical simulation that incoiporates gravitation and adiabatic cooling, we also find departure from analytical prediction in case of simulations 
with pre-heating. The non-gravitational processes that include pre-heating is not included in our analytical calculations. Thus departure from 
theoretical predictions provides a particularly interesting approach to separate out the effects of non-gravitational process on tSZ statistics. 
Non-Gaussianity contribution from different baryonic components: We have studied the contribution to the tSZ effect from individual compo- 
nents such as the "cold" gas T < lO'^K, uncollapsed moderate overdense "warm" gas 1 + S < 100 and the total contribution from all components 
that include the intra-cluster gas within the halos. We find that in addition to the total tSZ maps, the individual tSZ maps constructed using 
T < lO'^K, 1 + S < 100 components too can be described using our approach. For simulations the with pre-heating we find a significant departure 
for the contribution from "warm" gas component 1 -|- J < 100. The departure is more pronounced when the contributions from virialised halos 
are included. However we also find a near-perfect match for the contribution from the cold component. 

Many previous studies have considered a halo model based approach for modeling of the gas distribution in collapsed, virialized objects. In 
this approach, the specific number density, and the radial profile of these halos are analyzed using a Press-Schecther formalism or its variants. 
However, detailed modelling of the complete PDF or the bias is possible in this approach only in an order by order manner. In our extended 
peiturbative approach or the lognormal analysis, we go beyond the order-by-order approach and construct the entire cumulant generating function 
^{z) and relate it to that of the underlying density distribution (^{z). This allows us to reproduce and predict the entire PDF of the tSZ distribution 
for a specific smoothing angular scale. The statistical picture that we have developed here is complementary to that based on the Press-Schechter 
formalism. 

We would also l i ke to point out that some previ ous studies have modeled the statistics of tSZ effect using the hierarchical ansatz 
jValageas & sii3l 19991 : PValageas, Schaeffer & However in these studies the contribution from individual halos were computed us- 

ing a virialization scheme and equilibrium profile for the pressure distribution within the halos. The results presented here are complementary 
to such prescriptions as we focus on the large-scale distribution of ionized gas. Instead of modeling the contribution from individual halos, we 
directly link the baryonic pressure fluctuation responsible for the large-scale tSZ effect in terms of the underlying mass distribution. 

We have ignored the presence of instrumental noise in our results. A beam-smoothed noise of known PDF (Gaussian or otherwise) can always 
be convolved with the theoretical PDFs presented here before comparing them with any observational data. The method we pursue here also relies 
on having access to frequency cleaned tSZ maps. The tSZ effect can also be studied using cross-correlation techniques that involve external 
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tracers; such methods typically employ the "mixed" bispectrum. The results, however, lack frequency information and are typically dominated by 
confusion noise. Using frequency cleaned maps is expected to enhance the signal-to-noise significantly by exploiting frequency information, in 
the absence of which the background CMB plays the role of intrinsic noise that degrades the signal-to-noise ratio. It is also interesting to note that 
the removal of tSZ signals from the CMB maps may actually help detection of other sub-dominant effects. The study of PDFs presented here can 
play important role in this direction. 
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Figure Al. The left pane! depicts PDF p{5) as a function 1 + (5 for various values of ct = 0.1, 0.25, 0.5, 1.0, 1.5. Two different approximations are considered, the 
lognormal (long dashed) and the hierarchical ansatz (short dashed). The broader PDFs correspond to higher values of a. The bias b(> 5) is plotted as a function l + <5 
in the right panel. Two different approximations are considered, the lognormal distribution (long dashed) and the hierarchical ansatz (short dashed). The analytical 
results correspond to the large separation limit. 
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APPENDIX A: THE LOGNORMAL DISTRIBUTION 



The evolution of PDF of t he density field 5 can also be modelled using lognormal distribution. ( lHamiltonll985l : IColes & Jonesll99ll Bouchet et all 
1993: Kof mai l et aljl994l) . Detailed discussion of various issues involving of lognormal distribution can be found in'Sernardea u & KofmanHl995l) ~ 
Colombia 994^ We will use the following expressions for the PDFp((5) and the joint-PDF ^2) for our study (iTaruva et al .Il2002h : 



p{5)d5 = ^ 
E =ln(l + a^); 
p{Si, 52)dSidS2 = 



exp 



'2E2 



dS 
T+S' 



A = ln[(l + 5)\/(l + a2)] 
1 



^12 



exp 



S(A? 



2X12A1A2 



2(E2 



dSi dS2 

1 + (5l 1 + (52 ' 



A, = ln[(l + 5i)\/(lT^l; ^12 =ln(l+Cr(n,r2)). 



In the limiting case of large separation X12 — >■ we can write down the two point PDF as: 

p(Si,S2) = p(5i)p{S2)[l + b{5i)^f\ri,r2)b{S2)]; h{5. 



A./E. 



(Al) 
(A2) 
(A3) 

(A4) 
(A5) 



It is however easier to estimate the cumulative or integrated bias associated with objects beyond a certain density threshold 5q. This is de- 
fined as h{5 > So) — p{S)b{5)dS/ p{S)dS. In the low variance limit ct^ — > the usual Gaussian result is restored b{6) = 
The parameters A, Ai,A'i2,E that we have introduced above can be expressed in terms of the two-point (non-linear) correlation function 
^^^^(ri,r2) = ((5(ri)(5(r2)) = {Si52) and the nonlinear varia nce = (5^) of the smoothed de nsity field. Lognormal distribution has already 
been used to model the stati s tics o f weak lensing observables jMunshilbood : Irkruva et al.ll2002h and the clu stering of Lyman alpha absorption 
systems (e.g.iBi & Davidsoi] i ll 9971 )). and is known to model gravitational clustering in the quasilinear regime (iMunshi. Sahni. Starobinskvll 19941 : 
iMatarresse et al. 1992). In Figure JAU we compare the PDF and bias predictions from the lognormal model and hierarchical ansatz for various 
values of the variance. Though hierarchical ansatz and the lognormal model both predict nearly identical PDFs in the quasilinear regime, it is 
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p(<5) b((5) 

Figure Bl. A simple "flowchait" for a hierarchical approach is shown. It provides a summary of how the PDF p{5) and bias b{5) are constructed for a given model 
for the generating function for G{t). The construction involves the generating function for the normalised cumulants or Sp parameters and cumulant correlators 
Cpq i.e. (/'(z) and f^{z). These are related to p{S) and b{S) through Laplace transforms. The scaling functions h{x) and b{x) encodes the scaling properties of p{S) 
and b{S). We show the equation numbers that relate various quantities in the diagrams. The left hand side of the diagram correspond to one-point PDF and the right 
hand side correspond to the two-point PDF which is expressed through the bias function f'(i5). 



important to realize that the lognormal model is not a member of the family of hierarchical models, i.e. the PDF of lognormal model can't be cast 
into the form given in EQ.( l23t . 



APPENDIX B: HIERARCHICAL ANSATZ (MINIMAL TREE MODEL): A VERY BRIEF REVIEW 



As mentioned before the PDF p{S) and the bias b{S) can both be constructed from the knowledge of the VPF (j>{z) ~ X]p=i SpZ^ /p\ and its 
two-point analog t[z) = J]]^ Cpiz^ /p\. Where the parameters Sp and Cpi are normalized cumulants and cumulant correlators for the density 
field. 

The modelling of <j>[z) and t{z) needs a detailed knowledge of the entire coiTelation hierarchy. The detailed knowledge of the entire coiTela- 
tion hierarchy is encoded in the vertex generating function G{t)). Typically for large values of y the VPF exhibits a power law 4>{z) — az^^" . The 
parameter typically takes a value cj = .3 for CDM like spectra. There are no theoretical estimates of u and it is generally estimated from numerical 
simulations. For small but negative values of y the fu nctions ^(z) and t[z) develops a singularity in the complex plane which is described by the 
following parametrization dBalian & Schaeffenll989h : 



[z) = (f>s — asr{ujs){z — Zs) 



r{z) 



-bs{z- 



(Bl) 



The sing ularity structure of ^(z) and t{z) depends on the nature of the vertex generating function G(r) and its behaviour near the singularity Ts 
l lBalian & Schaeffer.l98A : 
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(B2) 



On the other hand the para meters lo and a can be des cribed in terms of a parameter Zs which in turn describes the exponential decay of the PDF 
for large density contrast 5 l lBalian & Schaeffeilll989h : 



_ kg -\-2 J ka/ka + 2 



1 1 {ka + 2y 

= Xi, — 



(B3) 



The PDF and the bias thus has two distinct regimes that are dictated by the two asymptotes. For inter mediate values of S the PDF shows a power 
law behaviour. The PDF and the bias are given by the following expression dBalian & Schaeffe3ll989h : 



p{5) 



1/2 



r[|(l + a;)] V 



1 + S 



{l-")/2 



For large values of 5 the PDF on the other hand shows an exponential behaviour dBalian & Schaeffeill 19891) : 
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1 + S 



exp 



1 + S 



1 + 



(B4) 



(B5) 
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At very small values of 5 the PDF shows an exponential decay which is described only by the parameter uj i lBalian&Schaeffeilll989l) : 



feW = -(^) (i^j ; v^{l + 5)[ar"''-"'\i2r"'-"''- (B7) 

The range of 5 for which the power law regime is valid depends on the value of ^2- For sm aller values of ^2 the power law regime is less 
pronounced. A more detailed discussion of these issues can be found in Munshi et alj 1 199%. The links to the gravitational dynamics in the 
quasilinear regime, for various approximations are discussed in lMunshi. Sahni. Starobinskvi l ll994l) . In this paper we have considered only a 
specific version of the hierarchical ansatz known as th e minimal tree model. Th i s is th e most popular version due to its simplicity; for variations 
and possible generalizations of minimal tree model see iBemardeau feSchaeffedd 19991) . 
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